R Notebook

library(tidyverse, verbose = F)
source("Funciones.R")

Lectura de Datos

Lo datos corresponde a una matriz 5x5x10. La metadata será X = 1…25 y Z = 1…10 (profundidad), ambos detectores tienen la misma metadata.

##Metadata
Metadata <- data.frame(x = rep(1:100, each = 11), 
                       z = rep(1:11, 100))

Comp.data <- data.frame(muestra = c("ss406", "ss407", "ss408", "ss409", "ss410"), 
                        C = c(0.19,0.50,0.28,0.11,0.39),
                        Mn = c(0.53,0.13,0.64,0.48,0.43),
                        Ni = c(1.69,0.61,4.58,3.14,2.04),
                        Cr = c(2.12,3.00,0.09,1.22,1.72),
                        Mo = c(1.03,0.82,0.14,0.77,0.41),
                        Cu = c(0.32,0.43,0.73,0.23,0.47))

DEMON

data.dir <- "Data/demon/"

wl <- read_tsv(paste(data.dir,"ss406.asc",sep = ""), col_names = F, progress = F, show_col_types = F) %>% 
        .[,1] %>% set_names("wl") %>% rowid_to_column()

## DEMON data
L_Demon <- map(c("ss406.asc","ss407.asc","ss408.asc","ss409.asc","ss410.asc"), 
               ~ read_tsv(paste(data.dir,.x,sep = ""), col_names = F, progress = F, show_col_types = F)) 

L_Demon <- L_Demon %>% set_names(c("ss406","ss407","ss408","ss409","ss410"))  

L_Demon <- L_Demon %>% map(~ .x %>% setNames(c("wl",paste("X", 1:(ncol(.x)-1), sep = ""))))

## elimina shot 1
L_Demon <- L_Demon %>% 
  map(~ .x %>% .[, c(1, which(Metadata$z != 1) + 1)]) ## elimina disparo #1 (limpieza)

g <- L_Demon %>% map(~ data.frame(.x[,1], Int = apply(.x[,2:ncol(.x)], 1, mean))) %>% bind_rows(.id = "id") %>% 
  ggplot(aes(x = wl, y = Int, color = id)) + geom_line() + labs(x = "Wavelength")
g %>% plotly::ggplotly()

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

## Normalizar
izq <- 267.70
der <- 267.75
int <- L_Demon %>% map_dfr(fun.int, izq = izq, der = der)
factor <-  1
int_norm <- int/factor

## acumular
int_norm <- sumar_filas_por_grupos(int_norm, n = 10, wl = F)

## Promediar
int_norm <- cbind(Comp.data, int_norm)
int_norm %>% pivot_longer(X1:X100, values_to = "Libs.Int") %>%
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean <- apply(int_norm %>% select(X1:X100), 1, promediar_grupos_aleatorios, n=10) %>% t() %>% data.frame()
int_mean <- cbind(Comp.data, int_mean)
int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_boxplot() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

MECHELLE

data.dir <- "Data/mechelle/"
## Mechelle data
L_Mechelle <- map(c("ss406.asc","ss407.asc","ss408.asc","ss409.asc","ss410.asc"), 
                  ~ read_tsv(paste(data.dir,.x,sep = ""), col_names = F, progress = F, show_col_types = F)) 

L_Mechelle <- L_Mechelle %>% map(~ .x %>% setNames(c("wl",paste("X", 1:(ncol(.x)-1), sep = ""))))
L_Mechelle <- L_Mechelle %>% set_names(c("ss406","ss407","ss408","ss409","ss410"))  

Exploracion Manual (Cromo)

## elimina shot 1
L_Mechelle <- L_Mechelle %>% 
  map(~ .x %>% .[, c(1, which(Metadata$z != 1) + 1)]) ## elimina disparo #1 (limpieza)

g <- L_Mechelle %>% map(~ data.frame(.x[,1], Int = apply(.x[,2:ncol(.x)], 1, mean))) %>% bind_rows(.id = "id") %>% 
  ggplot(aes(x = wl, y = Int, color = id)) + geom_line() + labs(x = "Wavelength")
g %>% plotly::ggplotly()

El spectrometer tiene un gran problema de calibracion. Los espectros se desplazan en longitud de onda. En el siguiente apartado implemento una funcion para recalibrar el espectro tomando como referencia una linea de la matriz (en este caso Fe 438.32).

Wavelength Re-Calibration

ref <- 438.3275
p <- data.frame(L_Mechelle[[1]][,1], Int = apply(L_Mechelle[[1]][,2:ncol(L_Mechelle[[1]])], 1, mean))
p$wl <- round(p$wl,4)
p <- which(p$wl == 438.3275)

L_Mechelle_new <- L_Mechelle %>% 
    map(~ fun.WL.calibration(old = .x, p = p))

L_Mechelle_new <- L_Mechelle_new %>% 
                    map(~ .x %>% drop_na() %>% filter(between(rowid, 1, 26818)) %>% select(!rowid))

g <- L_Mechelle_new %>% purrr::map(~ data.frame(wl = .x[,1], Int = apply(.x[,2:ncol(.x)], 1, mean))) %>% 
            bind_rows(.id = "id") %>% 
       ggplot(aes(x = wl, y = Int, color = id)) + geom_line() + labs(x = "Wavelength")
 
g %>% plotly::ggplotly()

Normalizacion: Estandar interno

¿Como tratar los datos? - ¿Que es una observacion?

• Metodo 1: Acumular (z) - Normalizar - Promediar (x)

## acumular
L <- L_Mechelle_new %>% map(~ .x %>% sumar_filas_por_grupos(n = 10))
  
## Normalizar
int <- L %>% map_dfr(fun.int, izq = 267.64, der = 267.74)       ## Cromo
factor <- L %>% map_dfr(fun.int, izq = 268.40, der = 268.50)    ## Fe
int_norm <- int/factor

## Promediar
int_norm <- cbind(Comp.data, int_norm)

int_mean <- apply(int_norm %>% select(X1:X100), 1, promediar_grupos_aleatorios, n=10) %>% t() %>% data.frame()
int_mean <- cbind(Comp.data, int_mean)

## graficar
int_norm %>% pivot_longer(X1:X100, values_to = "Libs.Int") %>% 
    ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
        geom_point()

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_violin() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

• Metodo 2: Normalizar - Acumular (z) - promediar (x)

## Normalizar: por espectro
int <- L_Mechelle_new %>% map_dfr(fun.int, izq = 267.64, der = 267.74)
factor <- L_Mechelle_new %>% map_dfr(fun.int, izq = 268.40, der = 268.50)
int_norm <- int/factor

## acumular: en Z
int_norm <- sumar_filas_por_grupos(int_norm, n = 10, wl = F)

## Promediar: en X
int_mean <- apply(int_norm, 1, promediar_grupos_aleatorios, n=10) %>% t() %>% data.frame()
int_mean <- cbind(Comp.data, int_mean)
int_norm <- cbind(Comp.data, int_norm)

## graficar
int_norm %>% pivot_longer(X1:X100, values_to = "Libs.Int") %>% 
    ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
        geom_point()

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_violin() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

Nota: En adelante una observacion sera el resultado de aplicar Metodo 1.

Normalizacion: suma total

Acumular (z) - Normalizar - Promediar (x)

## acumular
L <- L_Mechelle_new %>% map(~ .x %>% sumar_filas_por_grupos(n = 10))
  
## Normalizar
int <- L %>% map_dfr(fun.int, izq = 267.64, der = 267.74)       ## Cromo
factor <- L %>% map_dfr(fun.int, izq = 230, der = 1000)    ## Fe
int_norm <- int/factor

## Promediar
int_norm <- cbind(Comp.data, int_norm)

int_mean <- apply(int_norm %>% select(X1:X100), 1, promediar_grupos_aleatorios, n=10) %>% t() %>% data.frame()
int_mean <- cbind(Comp.data, int_mean)

## graficar
int_norm %>% pivot_longer(X1:X100, values_to = "Libs.Int") %>% 
    ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
        geom_point()

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_violin() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

Normalizacion: Porcion del espectro

## acumular
L <- L_Mechelle_new %>% map(~ .x %>% sumar_filas_por_grupos(n = 10))
  
## Normalizar
int <- L %>% map_dfr(fun.int, izq = 267.64, der = 267.74)       ## Cromo
factor <- L %>% map_dfr(fun.int, izq = 230, der = 320)    ## Fe
int_norm <- int/factor

## Promediar
int_norm <- cbind(Comp.data, int_norm)

int_mean <- apply(int_norm %>% select(X1:X100), 1, promediar_grupos_aleatorios, n=10) %>% t() %>% data.frame()
int_mean <- cbind(Comp.data, int_mean)

## graficar
int_norm %>% pivot_longer(X1:X100, values_to = "Libs.Int") %>% 
    ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
        geom_point()

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_point() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

int_mean %>% pivot_longer(X1:X10, values_to = "Libs.Int") %>% 
  ggplot(aes(x = Libs.Int, y = Cr, color = muestra)) + 
    geom_violin() + 
    geom_hline(yintercept = int_norm$Cr, lty = 2, col = "gray")

PLS (+Lasso) - Exploracion

library(mixOmics)